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Abstract 

Shapes of planar lipid monolayer domains at the air- water interface are theoretically and numer- 
ically investigated by minimizing the formation energy of the domains which consist of the surface 
energy, line tension energy, and dipole electrostatic energy. The shape equation which describes 
boundary curves of the domains at equilibrium state is derived from the first order variation of the 
formation energy. A relaxation method is proposed to find the numerical solutions of the shape 
equation. The theoretical and numerical results are in good agreement with previous experimen- 
tal observation. Some new shapes not observed in previous experiments are also obtained, which 
awaits experimental confirmation in the future. 
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I. INTRODUCTION 



In the last three decades lipid monolayer domains (LMDs) at an air-water interface have 
drawn much attention from biochemists and physicists. Experimental researchers have ob- 



served that the lipid monolayer domains dis 



S form, dumbbel 
with spiral arms 



n 



form, serpentine form[l|, 
2, 



ay many kinds of shapes, such as circular form, 
torus formj^, sl, fan form, multi-leaves form 



(also called labyrinthine form) and so on, which brings 
a theoretical question: How to understand these shapes? McConnell et al. proposed that 
the shape of LMDs was a natural result of competition between the energy due to the line 
tension of domain boundary and the dipole electrostatic energy of lipid molecules 
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By assuming that all lipid molecules align along the normal direction of the domain plane, 
they also proved that the total dipole electrostatic energy can be simplified as a double curve 
integral form along the domain boundary which will be called as McConnell energy in the 
following contents. Following McConnell's seminal work, Iwamoto and Ou-Yang proposed 
that besides the above two kinds of energy, the growth of a lipid domain should cost an ad- 
ditional surface energy due to the difference in the Gibbs free energy density between outer 
(hquid-like disorder) and inner (sohd-like order) phases |l2l]. Each observed LMD should 
minimize its " formation energy " (i.e., the sum of these three kinds of enrgy) 1^. Based on 
this formation energy, Iwamoto et al. explained the existence of several kinds of non-circular 



domains observed in the experiments 



HQ,Q. 



In fact, an implicit assumption to derive the McConnell energy is that the interaction 
between two dipoles is proportional to where ^ is the distance between them. It is this 
assumption that leads to the divergence of McConnell energy. As is well known, this as- 
sumption does not hold when £ is close to the dipole length(i.e., the length of lipid molecules 
or the thickness of LMDs). It is safe to use McConnell energy if the least distance between 
dipoles in the LMDs is much smaller than the thickness of LMDs. In the contrary case, the 
dipole electrostatic energy should be dealt with in the other way. Fortunately, when Langer 
et al. investigated the pattern formation in magnetic fluidjl^, as a byproduct, they proved 
that the total dipole electrostatic energy can also be simplified as a double integral form 
along the domain boundary which will be called as Langer energy in the following contents. 
The only difference between it and the McConnell energy is in the integrand. As a result, 
the Langer energy avoids the divergence of dipole electrostatic energy. In experiments, the 
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thickness of LMDs is about 2 nanometers while the least distance between dipoles in the 
LMDs might be several angstroms. Therefore, it is of significance to reconsider the shapes 
of LMDs based on the Langer energy. In this paper, the formation energy of a LMD is 
defined as the sum of the surface energy, the energy due to the line tension, and the Langer 
energy. The shape equation describing the boundary curve of the LMD is derived from 
the variation of the formation energy. We obtain many analytical and numerical solutions 
to the shape equation which agree well with the previous experimental results. We also 
present several kinds of shapes which have not been observed in experiments, which awaits 
the further conformation in the further experiments. 

The rest of this paper is organized as follows: In Sec. II, we present the formation 
energy of a LMD based on Langer energy, and then derive the shape equation of LMD by 
minimizing the formation energy. In Sec. Ill, we describe briefly our algorithm to solve 
the shape equation. In Sec. IV, we present many analytical and numerical solutions to the 
shape equation, and then compare them with experimental results. The last section is a 
brief summary and prospect. 



II. THEORETICAL MODEL 



A LMD with thickness h is schematically depicted in Fig. [T] where 2;— axis is the normal 
direction of the LMD and the projection of LMD in xy— plane is enclosed by a boundary 
curve {r(s)} with arc-length parameter s. t(s) and n(s) represent the tangent and normal 
vectors of point r(s) in the boundary curve, respectively. They satisfy dt/ds = nn and 
dn/ds = —K,t, where k is the curvature of the boundary curve. For simplicity, we introduce 
notations R^^ = r(/) — r(s), R = |Ris|, and R/s = His/R. The Langer energy (i.e., the dipole 
electrostatic energy) can be expressed as 1^: 

F,^J^f\i,f\i.I,l,I,M^lM!^, (1) 



2/^2 Jo Jo Jl Jl .fW+TI 



z 



/^2 



where / represents the integral along the boundary curve, is the dipole density in the 
LMD. 

The formation energy of the LMD can be expressed as 

F = APA + -f j^ds + FL, (2) 
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FIG. 1: Schematic illustration of the slab geometry of a lipid domain, where h, r(s), t(s), n(s) 
are respectively the thickness of the domain, the displacement vector of any point in the boundary 
curve, the tangent and normal vectors of the point in the boundary curve. 



where AP is the surface energy density introduced by Iwamoto and Ou-Yang in Ref. 12|. A 
is the surface area of LMD and 7 is hne tension of the boundary. 

Now we will drive the shape equation by using the variational method proposed in 
Ref. |l7i]. The variation in the tangent direction of the boundary curve will give a triv- 
ial identity. If we take the variation in the normal direction as 6r{s) = Q{s)n{s), it is not 
hard to obtain the following variational formulas: 

SA = ^rfsfi(s), (3) 



6 j>{ds) = - f dsK{s)n{s), (4) 



where " x " is the 2-dimensional cross product which satisfies a x b = 0162 — a2&i for 
any 2-dimensional vectors a = (01,02) and b = (61,62)- By considering Eqs.(l3])-(l5l), and 
lo dz dz'[R^ + {z - z'YY^/"^ = R~'^{y/WT~h^ - R),we can derive the shape equation 

AP - ^k{s) - (2/iV/i2) ^ dl[Ru X t(0] [\fl + {h/Kf - 1] = (6) 

from 5F = 0. 
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III. NUMERICAL METHOD 



In the last section, we have derive the shape equation ([6]) of LMDs by using variational 
method. However, it is very hard to obtain its solutions. Here we will develop a numerical 
method to solve it. In terms of the physical meaning of variational method, we can express 
the force in the boundary curve in a vector form as 

f(s) = - |aP - ^k{s) - (2/iV/i2) <fdl[Ris X t{l)][^/l + {h/Rf - 1]| n. (7) 

Thus Eq.([6]) reflects the force balance (f = 0) in the boundary curve. 

We introduce a virtual damping system evolving with dynamical equation 

^ = f(^). («) 

where r is virtual time. If a LMD whose boundary curve satisfies Eq.(l6]), its shape will keep 
unchange. On the contrary, f 7^ 0, the boundary curve will evolve according to Eq.®. In a 
long enough time, a damping system can usually approach to a state satisfying dr/dr = 
which implies f = 0. Thus the final configuration is naturally a solution to the shape 
equation ([6]). The above scheme to find numerical solutions of Eq.(l6]) is called the relaxation 
method, which had also been used to resolve hydrodynamics of the ferrofiuid drop pattern 

n n □ □ 

formation [16|, 118|, 119|, shape relaxation in L ang muir layers[20!] and hydrodynamics of the 



monolayer problem at the air-water interface 2l|] by Goldstein and his collaborators. 



In order to perform numerical calculations, we should transform the continuous equation 
([8]) into the discrete form by dividing the boundary curve into N segments with points 
and assuming the step of evolution time is tq. Then we label points in the boundary as 
1,2, - ■ ■ ,k — l,k,k + l, - ■ ■ ,N and the time sequence as tq, 2tq, ■ ■ ■ , (j — l)ro, Jtq, (j + 1)tq, ■ ■ ■. 
Eq.(IHl) is transformed into 



k „fc ck 



r„-Li - r 



j+i '-j - '-j 



ffro, (9) 



where is the discrete form of Eq.([7]). The geometric quantities in this equation can be 
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discreted as 



[ds]'; = |rf '-r^l, (10) 




(12) 



where can be obtained by rotating 



in 90° clockwise. 



The stop condition of our calculation of Eq.Q is Y.k=i kj+i < e for a large enough 
integer J and a small enough number e. 

IV. RESULTS 

A. ANALYTICAL SOLUTIONS 

1. Circular solutions 

For a circle with radius Rq, k = Because R^g and t(Z) are unit vectors, we can re- 

spectively express them as (cos 9is, sin 9is, 0) and (cos 9i, sin 6^^, 0). Eq. ([6]) is then transformed 
into: 



where x = a = , b = and E{.) is complete elliptic integrals of the second kind. 

From Fig [2], we find that a critical a* exists such that there are two solutions when 
< a < a* and no solution when a > a*. Additionally, there is only one solution for a < 
or a = a*. 

In the limit case of a; — > oo, we have 



ax + b = xE{2ix) — 2x^ 



(13) 



xE{2ix) - ^ - ln(8yex) 



(14) 



from which we easily obtain the radius of the circular domain 



7 



i?0 



8v^ 



(15) 



for vanishing AP. 



This result is similar to McConnell's result 
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6 



(16) 



with a cutoff 5 |lO|. 



6 




FIG. 2: Schematic diagrams of functions y = ax + h and y = xE{2ix) — 2x^. A critical a* exists so 
that there are two solutions when < a < a* and no solution when a > a*. Additionally, there is 
only one solution for a < or a = a*. 

2. Toroidal solutions 

Next, we investigate whether toroidal shapes are permitted in our model. The torus 



solution is mentioned in many references 
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13|. Here we only discuss the case of 



h <^ min{Ri — R2, R2}, where Ri and R2 are the outer and inner radii, respectively. 
From Eq. ([6]) we find the radii Ri and R2 should satisfy 



and 



A T^T^ 9i 8\/eRi n /■^'^ x(x — cos 6)d6 „ 

h Jo (x2 + 1 - 2xcos6')2 

2, 8^/eR2 2/"^^ a;(l -xcos^)^^ 



APR2 - 7 + uMn ^, - / ^ i ^dl = 0, (18) 

h Jo (x2 + l-2xcosr- ^ ^ 



2 



where x = R2/Ri- 

By adding Eq. (HZD to Eq. ([I8]), we obtain 

a{l + x) + \nx-2[E{k)-K{k)]=Q (19) 

where a = ^^t^, k'^ = — 4a;/(l — x)"^. K{k) and E{k) are respectively complete elliptic 
integrals of the first kind and the second kind, respectively. As shown in Fig. [3l there is 
a critical solution x ~ 0.182215 when a ~ 2.525115. This solution is close to the ratio 0.2 
observed in the experiment jsl. For a > 2.525115 there are two toroidal solutions. 



°5nf(x) 




FIG. 3: Schematic diagram of function /(x) = a{l + x) +lnx — 2[E{k) — K(k)] when a ~ 2.525115, 
there is a solution x = R2/R1 ~ 0.182215. 

Therefore toroidal shapes are indeed permitted under certain conditions in our modeh 
B. NUMERICAL SOLUTIONS 

There are four parameters, fi'^ , j , h and AP in our modeh We take the experimental 



values h ^ 2nm, /i^ = 5 x lO^^dyn, and 7 = 1.6 x IQ^'^dyn 22|. The surface energy density 
AP is an adjustable parameter, which can be understood as the Gibbs free energy difference 
between the fluid phase and solid phase [l^ . 

To make sure of our numerical method, we consider the circular domains. The radius of 
circular domains is about 20 fim in the experiment |22i]. Using the theoretical result Eq. (fUl) . 
we obtain AP ^ —5.05255 x 10~^dyn/cm. Adopting this parameter in our simulations, we 
obtain stable circular domain with radius about 20/im. 

Now we use our numerical codes to search for various shapes of domains evolving from dif- 
ferent initial configurations and AP, and then compare them with the experimental results. 
After the long time calculations, we obtain heart form (FigJlj), S form (Figl5]), moon form 
(FiglH]) and dumbbell form (Fig|7j) with adjustable parameters AP ^ 1.19338 x lO'^^dyn/cm, 
3.87582 X 10~'^dyn/cm, —1.42196 x I0~^dyn/cm, and 2.01274 x lO^^dyn/crru respectively. 



These shapes and their sized are in good agreement with the experiments [l, 7]. 

In addition, our model can also obtain one-lobed domain (also called labyrinthine pattern) 
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FIG. 4: Smooth heart. The top figure is cut from Fig. 5 (B) in Ref.[l|, where the bottom one is the 
numerical result with adjustable parameter AP ~ 1.19338 x 10~'^dyn/cm. 




FIG. 5: S domain. The figure in the top right hand corner is cut from the second picture of 
Fig. 3 in Ref.l?!, where the bottom left one is the numerical result with adjustable parameter 
AP w 3.87582 x W~^dyn/cm. 

and three-leaves domain shown in Figs. [S] and M under the parameters AP ^ —3.97255 x 
lO^'^dyn/cm and —1.66338 x 10~'^dyn/c'm, respectively. Their shapes and sizes are similar 
to those observed in Peter Kriiger and Mathias Losche's experiment l9|. 

n 

Furthermore, we not only obtain a six- leaves form, which occurs in the experiment [8|] in 
recent years, shown in Fig{TO]with adjustable parameter AP ^ —3. 11274 x 10~'^dyn/cm, but 
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FIG. 6: Moon. The right figure is cut from Fig. 5 (D) in Ref.[l|, where the left one is the numerical 
result with adjustable parameter AP « —1.42196 x lQ~^dyn/ cm. 
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FIG. 7: Dumbbell. The top figure is cut from Fig.l in Ref.[l|], where the bottom one is the 
numerical result with adjustable parameter AP ~ 2.01274 x 10~'^dyn/cm. 



also acquire some shapes with cusps observed in the experiments shown in Figs{TT] and [12] 
with adjustable parameters AP ^ 9.94745 x 10~^dyn/cm and AP ^ 9.72545 x 10~^dyn/cm, 
respectively. These results agree well with the experiments [l, 7]. 

Besides the shapes shown above, we can also produce some new stable shapes such as 
four-leaves, five-leaves, and six-leaves with cusps shown in FigJT3| which need to be verified 
by further biochemical experiments. 
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FIG. 8: One-lobed domain. The top figure is cut from Fig. 4 in Ref.[9|l, where the bottom one is 
the numerical result with adjustable parameter AP « —3.97255 x 10~^dyn/cm. 
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FIG. 9: Three- leaves form without spiral arms. The figure in the top right hand corner is cut 
from the fifth picture of Fig. 3 in Ref.[9], where the bottom left one is the numerical result with 
adjustable parameter AP « —1.66338 x 10~'^dyn/cm. 



The amazing aspect in our results is that we can obtain some stable shapes with "cuspidal 
points". Here the cusps are apparent ones viewed in the scale of micrometers. In our 
computational procedure, each cusps locally contain many discrete points in the scale of 
tens of nanometers such that the cusps are in fact the smooth curves. We find the average 
free energy per unit area of domains with cusps is less than those without cusps, that is. 
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FIG. 10: Smooth six-leaves form without spiral arms. The figure in the bottom right hand corner 
is cut from the second picture of Fig. 3 in Ref.[8], where the top left one is the numerical result 
with adjustable parameter AP « —3.11274 x 10~^dyn/cm. 
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FIG. 11: Bean with two cusps. The top figure in the bottom right hand corner is cut from the first 
picture of Fig. 3 in Ref.[7], where the top left one is the numerical result with adjustable parameter 
AP w 9.94745 x W-^dyn/cm. 

the former is stable than the latter. We show the bean shape, heart shape and multi-leaves 
shapes with cuspidal points in FigsiH], [12] and the last three pictures in Fig{T3l In the latest 
theoretical work [l^ , Iwamoto et al. analytically prove that the existence of the cuspidal 
points is reasonable and there are two kinds of cusps, one's curvature is positive, the other's 
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FIG. 12: Heart with a cusp. The figure in the bottom right hand corner is cut from Fig. 5 (D) in 
Ref.[l|, where the top left one is the numerical result with adjustable parameter AP ~ 9.72545 x 



FIG. 13: Multi-leaves forms. Top left, four-leaves form without cusps when AP ~ 1.63377 x 
10^ "^dyn/ cm. Top right, four-leaves form with cusps when AP « -2.06338 x 10""^ dyn/cm. Bottom 
left, five-leaves form with cusps when AP ~ —8.63377 x W~^dyn/cm. Bottom right, six-leaves 
form with cusps when AP ~ —1.34127 x W^^dyn/cm. 

is negative. However, their method is a local theory, in which one equilibrium domain shape 
can produce only one type of cusp and the number is also one. Our theory is a global theory, 
because we directly calculate the double curve integrals and do not use any approximate 
expansion method. Through numerical simulation, we further find that two kinds of cusps 
can coexist in one equilibrium domain, as shown in FigJTT| in which there are two kinds 




10 ^dyn/cm. 
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FIG. 14: Phase diagram with various AP (in the Log coordinate). 



of cusps. In fact, the domains with cusps have clearly been observed in serval experiments 



,13, 



23 



2J]. Thus the shapes with cuspidal points might be widespread in nature. 
We have investigated various shapes of domains under different AP. Now we draw a 
phase diagram in Fig. [HI It is found that the domain shape will become complex from 
the original circular domain when the parameter AP changes from —5.05255 x lO^^dyn/cm 
towards both ends of the coordinate axis. It seems that the domains with cusps generally 
occur in the narrow ranges of parameter AP. Two shapes of domains can coexist in a 
certain interval, such as heart shape and dumbbell shape as shown in FigsJT^ which agrees 

that reveals some shapes can coexist under some 

parameters. 

V. CONCLUSION 

In the above discussions, we have theoretically and numerically investigated the shapes 
of planar lipid monolayer domains by minimizing the formation energy of the domains. 
We obtain some shapes, such as circular form, S form, dumbbell form, serpentine form 
and achiral multi-leaves form and so on. Most of them are observed in the experiments 

mm 

. The four-leaves, five-leaves, and six-leaves forms await the future experimental 
conformations. 

Our model might be extended to explain the shapes of rafts [25] consisting of sphingolipids 
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and cholesterol in the cell membranes. Although the cell membrane is a bilayer, the raft 
domains are generally in one of the monolayers. There are net dipoles in the raft domains, 
thus our model might be applied to describe the shape of raft domains. 
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